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Juring  tne  past  two  years  a  very  substantial  amount  of  progress  lias  been  made  toward  the 
main  goal  of  this  project,  which  is  to  calculate  the  activation  energy  of  explosive 
molecules.  The  progress  has  come  in  the  form  of  four  different  approaches  to  the  problem: 

i  Cl  method;  2)  Gaussian  82  Computer  Program;  3)  MNDOC  procedure;  4)  Green's  function 
techniques.  We  have  recently  obtained  Dr.  Henry  F.  Schaefer's  Cl  programs,  which  are 
generally  considered  to  be  the  best  and  fastest  Cl  programs  in  the  world.  We  have  almost 
finished  adapting  these  programs  to  the  Cray  computers,  which  arc  among  the  fastest 
umpjtcrs  in  the  world  and  are  often  referred  to  as  super  computers.  When  this  adaptation 
is  completed,  we  will  have  the.  first  Cray  version  of  Dr.  Schaefer's  Cl  programs,  and  in 
uoditiun,  we  have  access  to  several  large  Cray  computers  on  which  to  run  the  programs. 
r..4  u  result,  we  expect  to  perform  very  large  and  highly  accurate  Li  calculations  on 
cAf.lu-.ive  molecules  of  interest.  The  new  Cray  version  of  these  Cl  programs  and  the 
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cotibidCMea  to  be  the  state  of  the  art  program  in  the  area  of  Mol  1 er-Pl esset  perturbation 
theory,  we  have  recently  obtained  access  to  a  new  Cray  version  of  Gaussian  82  which  is 
not  yet  available  to  the  general  scientific  community.  With  the  combination  of  this  new 
program  and  tne  large  cray  computers  to  which  we  have  access,  we  expect  to  perform 
activation  energy  calculations  that  are  better  and  more  accurate  than  those  that  were 
previously  possible  using  Gaussian  82.  The  MNDOC  method,  which  is  a  new  correlated 
version  of  the  MNLiQ  method,  has  been  used  to  compute  accurate  values  for  the  activation 
energy  of  methyl  nitrite.  We  have  recently  used  Green's  function  theory  to  derive  a  new 
one-electron  equation  that  goes  beyond  previous  one-electron  equations  to  include  higher 
order  correlation  terms.  As  correlation  plays  a  very  important  part  in  molecular  reactions, 
tnis  new  equation  is  expected  to  lead  to  significant  improvements  in  the  calculation  of 
.molecular  activation  energies. 
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Summary  of  Accuiiipl  i  ^huie/iLs 

During  the  past  two  years,  a  very  substantial  amount  of  progress  has  been 
made  toward  the  main  goal  of  this  project,  which  is  to  calculate  the 
activation  energy  of  explosive  molecules.  The  progress  has  come  in  the 
form  of  four  different  approaches  to  the  problem: 

1)  Cl  Calculations 

It  is  well  known  that  a  full  Cl  qives  the  best  possible  molecular 
total  energies  and  activation  energies  that  can  be  computed  for  a 
given  basis  set.  However,  since  Cl  calculations  can  be  very  time- 
consuming,  one  must  have  very  fast  programs  and  very  large  computers 
in  order  to  make  the  calculations  feasible  for  explosive  molecules 
of  interest.  We  have  recently  obtained  Dr.  Henry  F.  Schaefer's  Cl 
programs,  which  are  generally  considered  to  be  the  best  and  fastest 
in  the  world.  In  addition,  we  have  access  to  several  large  Cray 
computers,  which  are  among  the  fastest  in  the  world  and  are  often 
referred  to  as  super  computers.  We  have  almost  finished  converting 
Dr.  Schaefer's  Cl  programs  to  the  Cray.  Once  we  get  these  progrrms 
running  on  the  Cray,  we  expect  to  perform  very  accurate  Cl  calcula¬ 
tions  on  many  explosive  molecules  of  interest.  The  new  Cray  version 
of  these  Cl  programs  and  the  calculations  that  will  be  done  with  it 
are  expected  to  arouse  a  large  amount  of  interest  in  the  scientific 
community. 

2)  Gaussian  82 

The  Gaussian  82  computer  package  is  generally  considered  to  be  the 
"state  of  the  art"  program  in  the  area  of  Mol ler-Plesset  perturbation 
theory.  This  type  of  perturbation  theory,  if  carried  to  high  enough 
order,  can  yield  very  accurate  results  for  molecular  activation 
energies.  We  have  very  recently  obtained  access  to  a  new  Cray  version 
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of  Gaussian  82  which  is  not  yet  available  to  the  general  public. 

With  the  combination  of  this  new  proqram  and  the  large  Cray  computers 
to  which  we  have  access,  we  expect  to  perform  activation  energy 
calculations  that  are  better  and  more  accurate  than  those  that  were 
previously  possible  using  Gaussian  82. 

3)  MNDOC 

The  MNDOC  method  is  a  semi-empirical  technique  which  is  very  fast 
and  could  therefore  be  used  to  compute  the  activation  energies  of 
large  explosive  molecules.  One  indication  of  MNDOC's  potential  is 
that  it  has  yielded  accurate  values  for  the  activation  energy  of 
methyl  nitrite.  Further  improvements  are  needed  in  this  method 
before  it  can  be  applied  to  a  wide  range  of  explosive  molecules. 

We  have  recently  discovered  a  modification  that  may  greatly  improve 
the  method's  accuracy.  This  modification  consists  of  adding  a  two- 
configuration  SCF  to  the  MNDOC  programs.  The  two-configuration  SCF 
is  expected  to  work  well  for  molecules  in  which  the  ground  state 
and  first  excited  state  are  close  in  energy. 

4)  Green's  Function  Techniques 

We  have  recently  used  Green's  function  theory  to  derive  a  new  one- 
electron  equation  that  goes  beyond  previous  one-electron  equations 
to  include  higher  order  correlation  terms.  As  correlation  plays  a 
very  important  part  in  molecular  reactions,  this  new  equation  is 
expected  to  lead  to  significant  improvements  in  the  calculation  of 
molecular  activation  energies. 
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Introduction 

The  goal  of  this  project  is  to  calculate  the  amount,  of  shocking  or 
jostling  that  a  solid  explosive  can  withstand  before  detonating.  The  ultimate 
objective  is  to  determine,  by  means  of  the  calculations,  how  to  modify  the 

explosive  so  that  it  still  retains  all  of  its  explosive  power  but  becomes 

less  sensitive  to  shock.  These  less  sensitive  but  equally  powerful  explosives 
would  reduce  the  number  of  accidents  that  occur  among  personnel  involved  in 
handling,  storing,  transporting,  and  using  the  explosives. 

The  detonation  process  is  complicated  and  involves  many  factors, 
including  "hot  spots",  propagation  of  acoustic  or  detonation  waves,  inter¬ 
action  (or  scattering)  between  detonation  waves  and  explosive  molecules,  etc. 
However,  the  key  factor  in  an  explosive's  sensitivity  to  shock  is  the 
activation  energy  of  the  molecules  that  make  up  the  explosive.  When  an 
explosive  receives  a  shock,  the  energy  of  the  molecules  contained  in  the 

explosive  is  increased.  In  the  region  of  the  shock,  if  the  energy  of  a 

sufficient  number  of  molecules  is  raised  above  the  activation  energy,  then 
the  detonation  begins.  Thus  if  one  can  increase  the  activation  energy,  a 
bigger  shock  will  be  required  for  detonation  and  the  explosive  will  be  safer 
to  handle,  "herefcre  the  first  step  in  the  study  of  detonations  is  to 
determine  accurate  values  for  the  activation  energy  of  the  explosive  molecules. 

One  can  determine  the  activation  energy  of  an  explosive  molecule  by 
calculating  the  total  energy  of  the  molecule  as  a  function  of  a  reaction 
coordinate  such  as  a  bond  length  or  bond  angle.  At  the  value  of  the 
reaction  coordinate  where  the  total  energy  goes  through  a  maximum,  dissocia¬ 
tion,  or  detonation,  occur:.,  and  the  activation  energy  is  the  difference 
between  this  maximum  total  energy  and  the  total  enerqy  of  the  molecule  in 
its  equilibrium  state.  Thus  in  order  to  obtain  accurate  values  for 


molecular  activation  energies,  one  must  be  able  to  compute  accurate  values 
for  molecular  total  energies. 

It  is  well  known  that  in  order  to  obtain  accurate  and  useful  values 

for  molecular  total  energies,  it  is  essential  to  include  correlation  effects. 
.•Je  have  thoroughly  studied  four  different  methods  for  including  correlation  i 

molecular  calculations: 

1 )  Cl  Method 

2)  Gaussian  82  Computer  Program 

3)  MNDOC  Procedure 

4)  Green's  Function  Techniques 

Very  substantial  progress  has  been  made  in  all  four  of  these  areas  and  they 
all  show  great  promise  for  yielding  accurate  molecular  activation  energies. 
Eacti  of  these  methods  is  described  in  detail  below. 
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I.  Cl  Method 


As  mentioned  above,  a  very  attractive  feature  of  a  full  configuration 
interaction  (Cl)  calculation  is  that  it  will  yield  the  best  possible  molecular 
activation  energy  that  can  be  computed  for  a  given  basis  set.  Often  it  is 
not  necessary  to  carry  the  Cl  all  the  way  to  completion,  but  rather  to 
converge  the  Cl  to  whatever  degree  of  accuracy  one  desires  for  the  molecular 
activation  energies.  However,  for  explosive  molecules  of  interest,  an 
accurate  converged  Cl  can  be  very  time-consuming  and  therefore  one  must  have 
very  fast  computer  programs  and  very  fast  computers  in  order  to  make  the 
calculation  feasible. 

We  have  recently  obtained  Dr.  Henry  F.  Schaefer's  Cl  programs^  which 
are  generally  considered  to  be  the  best  and  fastest  in  the  world.  In 
addition  to  acquiring  these  programs,  we  are  at  present  making  a  major 
improvement  in  their  usefulness  by  adapting  them  to  the  Cray  computers,  which 
are  among  the  fastest  computers  in  the  world  and  are  often  referred  to  as 
super  computers.  When  this  adaptation  is  completed,  we  will  have  the  first 
Cray  version  of  Dr.  Schaefer's  Cl  programs,  ana  in  addition,  we  have  access 
to  several  large  Cray  computers  on  which  to  run  the  programs.  As  a  result, 
we  expect  to  perform  very  accurate  Cl  calculations  on  explosive  molecules  of 
interest.  The  new  Cray  version  of  these  Cl  programs  and  the  calculations 
that  will  be  done  with  it  are  expected  to  arouse  a  very  large  amount  of 
interest  in  the  scientific  community. 

In  addition  to  their  speed,  another  important  feature  of  Dr.  Schaefer's 
Cl  programs  is  their  ability  to  perform  a  two-configuration  SCF  and  then 
do  a  Cl  from  that.  While  the  starting  point  for  Cl  calculations  is  usually 
a  single-configuration  SCF,  or  Hartree-Fock  calculation,  sometimes  the  Cl  will 
converge  significantly  faster  if  the  starting  point  is  a  two-configuration 
SCF.  This  is  particularly  true  for  molecules  like  ni tromethane  in  which  the 


ground  state  and  first  excited  state  are  close  in  energy,  and  may  also  be 
true  for  larger  explosive  molecules  that  contain  the  nitro  group,  such  as 
nitrobenzene  and  trinitrotoluene  (TNT).  Thus  the  two-configuration  SCf 
further  increases  the  ability  of  Dr.  Schaefer's  programs  to  yield  a  converged 
Cl  for  explosive  molecules  of  interest. 

11  12 

We  now  give  a  brief  description  of  the  Cl  method  ’  .  The  Cl  procedure 

is  used  to  solve  the  time- i ndependent ,  nonrelativistic,  rigid-nuclei, 
electronic  Schroedinger  equation: 

Hy  =  Ey  (1) 

where  the  Hamiltonian  H,  in  atomic  units,  is  given  by 

H  =  1'  ,  +  7.  h  f  c  g  (?-) 
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Here  the  upper-case  Latin  subscripts  represent  the  (rigid)  nuclei,  the  lower¬ 
case  Greek  subscripts  represent  the  electrons,  the  Z’s  are  atomic  numbers, 
and  the  r's  (t'AR,  rft,  ,  r  )  are  the  distances  between  the  particles.  The  VNfJ 
term  is  a  constant  since  the  nuclei  remain  fixed,  and  therefore  we  can  ignore 
it  throughout  the  calculation. 

In  order'  to  solve  the  Schroedinger  equation,  we  write  tire  wave  function 
,  as  a  linear  combination  of  symmetry  adapted  configuration  functions  : 


where  the  1 


s  are  linear  combinations  of  Slater  determinants.  In  a  Cl 
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calculation,  one  need  include  only  conf igurations  of  the  same  symmetry  as  the 
particular  electronic  state  being  investigated  (the  reference  state),  since 
Hamiltonian  matrix  elements  between  configurations  of  different  symmetry  are 
zero.  Hence  the  configuration  functions  are  chosen  to  have  the  symmetry  of 
the  reference  state.  Since  the  total  number  of  configurations  of  all 
symmetries  can  be  very  large,  the  symmetry  adaptation  of  the  configuration 
functions  greatly  reduces  the  size  of  the  Hamiltonian  matrix  and  simplifies 
the  calculation. 

The  coefficients  c$  are  chosen  to  minimize  the  energy 

E  =  illHjli  (7) 

according  to  the  variation  principle.  The  variation  principle  ensures  that 
the  lowest  calculated  energy  is  an  upper  bound  for  the  exact  lowest  eigenvalue 
of  H  (i.e.,  the  exact  ground  state  energy).  It  also  ensures  that  the  lowest 
calculated  energy  for  each  symmetry  is  an  upper  bound  to  the  exact  lowest 
eigenvalue  of  that  particular  symmetry.  Thus  if  one  is  investigating  an 
excited  state  which  is  the  lowest  of  a  particular  symmetry  type,  one  need 
not  include  any  lower  states  in  the  same  calculation. 

Application  of  the  variation  principle  leads  to  the  well-known  eigenvalue 
e  :  j  a  t  i  o  n 

l  = 0  (8) 

w'ich  is  solved,  by  standard  techniques,  for  the  eigenvalues  E  and  coefficients 
c*.  In  the  above  equation,  H  and  are  the  Hamilton, an  and  overlap  matrix 
elements,  respectively,  and  are  given  by 
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(10) 
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TKe  main  computational  step  in  a  Cl  calculation  is  t Vie  evaluation  of  the 
Hamiltonian  matrix  elements  H^. 

In  order  to  construct  the  configuration  functions  and  evaluate  the  matrix 
elements,  one  chooses  a  set  of  atomic  basis  functions  .  One  then  writes  the 
atomic  orbitals  as  a  linear  combination  of  the  basis  functions 


A.  =  l  v  U  • 
V1  n  XP  PI 


01) 


where  the  coefficients  U  .  are  determined  by  an  SC.'  (self-consistent  field)  or 

pi 

similar  calculation.  The  evaluation  of  the  Hamiltonian  matrix  elements  is 
greatly  simplified  if  one  chooses  the  atomic  orbitals  <f.  to  be  orthonormal. 

T fie  Hamiltonian  matrix  elements  are  then  given  by 


1  ,  =  z  asth .  .  +  t  bs^,  g .  .. 

st  jj  'J  ’J  ijU  Uu9Uki 


(12) 


where 
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ar.  1  the  coefficients  as^  and  b^,  are  determined  by  a  projective  reduction 
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The  transformation  of  the  two-electron  integrals,  fig.  (16),  can  be  quite  time- 
consuming  and  must  therefore  be  handled  efficiently. 

To  summarize,  the  major  computational  steps  in  a  Cl  calculation  are  the 
fol lowing : 

1)  Choose  a  basis  set  and  compute  the  basis  set  integrals 

2)  Determine  the  orbitals  by  means  of  an  SCF  calculation 
2)  Transform  the  basis  set  integrals  to  orbital  integrals 

4)  Choose  and  construct  a  set  of  symmetry  adapted  configuration 
functions  for  the  electronic  state  that  is  being  investigated 

5)  Calculate  the  Hamiltonian  matrix  elements  with  respect  to  the 
configuration  functions 

6)  Compute  the  lowest  eigenvalue(s)  and  eigenvec tor( s )  of  the 
Hamil ton i an  matrix 

As  mentioned  previously,  a  full  Cl  calculation  (one  in  which  all  the  configura¬ 
tions  of  a  particular  symmetry  are  included)  yields  the  best  possible  total 
energy  that  can  be  computed  for  a  given  set  of  basis  functi  as  y  . 


II.  Ga us s  i a n_  82  Compu ter  Prog  t  din 

The  Gaussian  81’  computer  program  is  the  latest  in  a  series  of  moiecular 
ertital  programs  by  Joiin  Pople  and  his  group  at  Carneg i e-Mel  1  on  University. 
Gaussian  88  is  generally  considered  to  be  the  "state  of  the  art"  program  in 
the  area  of  perturbation  theory,  and  as  such  it  has  the  potential  t/  yield 
very  accurate  values  for  molecular  activation  energies.  As  in  the  case  of 
the  Cl  calculations  discussed  previously,  the  perturbation  theory  calculations 
can  be  very  time-consuming  and  therefore  very  fast  computers  are  needed  to 
obtain  accurate  results  for  explosive  molecules  of  interest.  We  have  recently 
made  substantial  progress  in  this  area  by  gaining  access  to  a  new  Cray 
.ersion  of  Gaussian  82  that  is  not  yet  available  to  the  general  public.  In 
audition,  we  have  access  to  several  large  Cray  computers  on  which  to  run  the 
program.  With  these  facilities  we  expect  to  calculate  activation  energies 
that  are  better  and  more  accurate  than  those  that  have  previously  been  computed 
w  •  *  h  thi  Gaussian  :>2  method. 

fne  of  the  main  features  of  Gaussian  82  is  its  capability  to  perform 
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, !  1  er-l-iossct  perturbation  theory  calculations  up  to  fourth  order.  This 

i  :  s  ' g n i f  i  r  ! n t  improvement  over  the  previous  version  of  the  program, 

••  "  ,  . , ; : i i, ! ,  could  carry  Mol  1  er-Plesset  perturbation  calculations  up 

*  .<■  •  !  e’er,  fur  many  molecules,  fourth  order  Moll  er-Plesset  perturbation 

ta  :  I  ,n  ge  part  of  the  correl  at  ion  energy  and  hence  yields  accurate 
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,>i;‘  ■;  energies.  In  addition,  Gaussian  82  can  perform  coupled  cluster 

1  ,i  P  i  ,  r;  ,  un  to  second  oi  ier.  The  coupled  cluster  method  also  lias  the 
<  ip.:  i  *  /  t  .  yield  accurate  activation  energies  for  explosive  molecules. 


B.  SECH  relation 

Carrying  out  stop  1)  of  the  iterations,  wo  put  ::(J  into  Eq.  (13).  Then 
Eq.  (13)  becomes  the  one-electron  Hartree  equation: 

[h( x ) +V(x ) ]Ug . (x )  =  E0iuQi(x)  (i; 

and  Gg  is  given,  according  to  Eq.  (12),  by 

,  .  ..  u0i(x')u0ilx> 

Grt(x,X,a>)  —  z.  r  .  ( 1  £ 

U  i  °~E0i 

tutting  ;.g  into  Eg.  (2),  we  obtain  the  correspond i ng  zeroth  order  vertex 


tunc t  i  o  n 


: nvl 2,3)  =  ‘(12)5(13)  . 


"hern  putting  ;n  and  GQ  into  Eq.  (3),  we  obtain  the  zeroth  order  polarization 


2)  =  -iG0(l?)fi  (21) 


wuicn,  L/  Eqs.  (4)  and  (3),  yields  the  zeroth  order  dielectric  function  and 


greened  interaction: 


(  1  2  !  -  i :>  f 1  \  ■  <  ( -to\  <  >.  ■ 

\  J  J  •  r\\  •  J  *  \  *-  /  >  \  t 


wf)(  12)  =  /v  (13)  cE  (  32 )d ( 3 )  . 


The  above  formula  for  the  polarization,  tq.  (20),  is  a  well-known  result 
uallcd  the  random  phase  approximation  (RBA)  !  'w,  and  :  and  Wg  are  the  RPA 
dielectric  function  and  screened  interaction: 


PP.?fS  1  <■ ) 


1  2 ) 


,..(i2)  ,  „;i 


Eqs.  ( 1 ) - ( 5 )  can  bo  iterated  to  obtain  successively  more  accurate 
expressions  for  Z,  !'  and  V  as  functionals  of  G.  At  cacti  stage  of  t tie 
iterations,  the  Green's  function  is  determined  by  z  through  Eqs.  (12)  and  (13). 
We  start  the  iterations  by  setting  z  equal  to  zero: 

Z0  =  0  .  (16) 

The  iterations  are  then  carried  out  in  the  following  sequence: 

1)  determines  Gg  through  Eqs.  (12)  and  (13) 

2)  Xg  and  Gg  determine  Tg  through  Eq.  (2) 

3)  Tg  and  Gg  determine  Pg  through  Eq.  (3) 

4)  Pg  determines  eg  and  Wg  through  Eqs.  (4)  and  (5) 

3)  ig  and  Wg  determine  X^  and  G^  through  Eqs.  ( 1  ) ,  ( 1 2 )  and  (13)  (i.e.,  1'g 
and  Wg  arc  substituted  into  Eq.  (1)  and  then  Lqs.  (1),(12)  and  (13)  are 
solved  sel f-consistently  for  Z  ^  and  G^ ) 

•'■)  :  i  and  G,  determine  :  ^  through  Eq.  (2) 

’)  . i  and  G|  determine  P-|  through  Eq.  (3) 

■3)  P,  determines  ,  and  W,  through  Eqs.  (4)  and  (3) 

i'i  and  Wj  determine  20  and  through  Eqs.  ( 1 ) ,  ( 1 2 )  and  (13)  (i.e.,  :  -j 
and  Wj  are  substituted  into  Eq.  (1)  and  then  Lqs.  (1),(12)  and  (13) 
sol  ved  sel  f-cons i stentl y  for  and  Go) 


a  re 
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0 )  =  (x1  ,t] ) 


(1  )  =  (x1 ,t1 +n) 


where  x  includes  space  and  spin  coordinates  (x)  =  (r,c)  and  is  understood 


The  Green's  function  satisfies  the  equation 


[  -  h(x1  )-V(l  )3G(12)  -  Jv(13)G(32)d(3)  =  6(12)  (1C 

where  h  is  the  kinetic  energy  plus  the  interaction  of  the  electron  with  the 
nuclei,  and  V  is  the  average  potential: 


V(l)  =  -i /v(13)G(33f)d(3) 


fhe  Fourier  transform  of  the  Green's  function  with  respect  to  time  is  given 


U.(x’)u.(x) 
G(x  ,x  '  ,u>)  =  L  - -f - 


where  is  the  frequency,  and  the  amplitudes  u.  and  energies  E-  satisfy  the 
following  equation  (which  is  the  Fourier  transform  of  Cq.  (10)): 


[h( x )  +V( x ) ]u  .  (x )  f  /:;(x,x  '  )ui  ,x  '  )dx  ’  i  .u  . 


V(x)  =  /  v  ( r , r ' )p(x')dx' 


1  occ . 


u  i  (  X  )  U  j  ( x )  . 


If  the  self-energy  y  is  independent  of  frequency,  then  the  amplitudes  and 
energies  given  by  Eq.  (13)  are  one-electron  wave  functions  and  energies,  and 
unc  ran  work  in  the  one-electron  picture. 
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A.  Green's  Function  Formalism 


For  a  neutral  system  containing  N  electrons,  one  can  use  Green's  function 

21-9"? 

theory  to  obtain  the  following  set  of  self-consistent  equations  ~v  for  the 
self-energy  I,  vertex  function  r,  polarization  P,  dielectric  function  e,  and 


screened  interaction  W: 

1(12)  =  i /W(1  +3 ) G(  1 4 ) I'(42,3)d(34  )  (1) 

r(l  2,3)  =  6(12)6(13)  +  G(46)r(67,3)G(75)d(4567  )  (2) 

P(  1  2 )  =  -i/G(23)r(34,l )G(42  )d(34 )  (3) 

c(l 2)  =  6(12)  -  / P ( 1 3)v(32)d(3)  (4) 

W(12)  =  /v(l 3)e"1 (32)d(3)  (5) 

w'  rTe  v  is  the  bare  instantaneous  Coulomb  interaction: 


v ( 1  2)  =  v (r-j  ,r2)6(t1  ,t2) 
2 

v(rl  ,r2)  = 


'  rl ”r2 


(6) 

(7) 


the  one-electron  Green's  function: 


G ( 1  2 )  =  -i[--;(l)y+(2)>0(trt2)  -  •-,"(;'),(;)  -(t;,-t.)] 


(8) 


:rt2) 


i  tr-t2 
0  t^  ■  t2 


(9) 


e  is  the  charge  of  an  electron,  if  is  the  Heisenberg  field  operator,  and  the 
brackets  •  •  indicate  the  ground  state  expectation  value.  In  the  above 

equations,  we  have  used  the  notation 
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above  is  the  screened  interaction  used  in  the  expansion  of  the  self-energy. 
This  screened  interaction  greatly  speeds  up  the  convergence  of  t lie  series 
for  the  self-energy. 

The  first  term  in  the  series  for  the  self-energy  leads  to  the  screened 
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exchange  plus  Coulomb  hole  (SECH)  one-electron  equation  which  is  solved 

sel f-cons i st entl y  to  obtain  correlated  molecular  energies.  The  SECH 

equation  has  never  been  used  before  in  molecular  calculations  and,  as 

mentioned  above,  is  expected  to  yield  accurate  values  for  molecular  energies. 

The  second  term  in  the  self-energy  series  leads  to  a  new  one-electron 

equation  that  goes  beyond  the  SECH  equation  to  include  higher  order  correla- 
27 

tion  effects.  This  new  equation,  wiiich  we  have  called  the  SECH2  equation, 
has  never  been  used  before  in  either  molecular  or  solid  state  calculations. 

It  is  a  completely  new  equation  that  contains  even  more  correlation  than  the 
SECH  equation  and  is  therefore  expected  to  produce  significant  improvements 
in  the  calculation  of  correlated  molecular  energies.  The  basic  Green's 
fit. (tion  formalism,  the  SECH  equation,  and  the  SECH2  equation  are  described 
in  detail  in  Sections  A-C. 


IV.  Green  '  s  J  une t  ion  Techniques 

v  1  _  2  j 

We  have  developed  a  new  formal  ism"  ,  based  on  Green's  function  theory, 
for  including  correlation  in  molecular  calculations.  To  our  knowledge,  this 
approach  has  never  been  applied  to  molecules.  In  this  new  formalism,  one 
uses  Green's  function  theory  to  obtain  a  set  of  successively  more  accurate 
expressions  for  the  self-energy,  which  contains  all  of  the  exchange  and 
correlation  effects.  Using  this  new  formalism,  one  can  include  correlation 
effects  in  a  successively  more  accurate  way,  until  one  obtains  the  desired 
degree  of  accuracy  for  the  molecular  energies.  In  this  process,  new 
molecular  one-electron  equations  will  be  developed. 

One  of  the  main  advantages  of  this  Green's  function  formalism  over 
other  perturbation  type  approaches  is  that  the  self-energy  is  expanded  in 
powers  of  a  screened  interaction  rather  than  a  bare  Coulomb  interaction.  The 
screening  of  the  interaction  greatly  speeds  up  the  convergence  of  the  power 
series  for  the  self-energy. 

Using  only  the  first  term  in  this  series  for  the  self-energy,  we 

2i  ?3 

obtained  excellent  results  for  solids."  The  question  is,  can  this 

Green's  function  method  also  be  applied  to  smaller  systems  such  as  molecules? 
We  now  present  very  strong  evidence  that  Green's  function  methods  can  be 
successfully  applied  to  molecules.  There  is  a  very  large  literature  on 
the  use  of  Green's  function  approaches  to  calculate  correlated  molecular 
energies.  References  28-85  refer  to  some  of  the  authors  who  have 
successfully  applied  Green's  function  techniques  to  molecules.  These 
papers,  which  are  only  a  partial  list  of  the  large  body  of  publications  in 
this  area,  clearly  show  that  Green's  function  approaches  can  be  successfully 
a ppl  i ed  to  mol ecul es . 

Once  again,  it  she  aid  be  noted  that  what  distinguishes  the  Green's 
function  method  of  this  proposal  from  the  Green's  function  methods  cited 


is  true,  then  including  more  excitations  wuuld  increase  the  error  produced  by 
this  inconsistency.  This  would  explain  why  BWLN2,  which  has  single,  double, 
triple  and  quadruple  excitations,  gives  worse  results  than  BWL'Nl  ,  which  has 
only  double  excitations. 

The  BWEN  method,  which  is  fast  and  can  thus  be  applied  to  relatively 
large  molecules,  yields  good  results  for  the  activation  energy  of  methyl 
nitrite.  The  results  are  not  as  good  for  nitromethane  and  thus  some 
improvements  are  needed  before  the  BWEN  procedure  can  be  generally  applied  to 
larger  molecules.  One  possible  improvement  is  the  incorporation  of  a  two- 
configuration  SCF  into  the  MNDO  method.  The  two-configuration  SCF,  rather 
than  a  single-determinant  SCF  (HF),  can  then  be  used  as  the  zero-order 
reference  state  and  the  starting  point  for  excitations  in  the  BWEN  second- 
u"der  perturbation  treatment.  The  two-configuration  SCF  is  expected  to  work 
we i 1  for  molecules  such  as  nitromethane,  in  which  the  ground  state  and  first 
excited  state  are  close  in  energy. 


Reac  tion 

Experiment 

MNDO/CI 

BWEN 

BH  L  Ml 

BWEf 

CH,N0„  -  CH, 

J  L  J 

+  no2 

59.0 

46.5 

48.2 

35.8 

29.1 

Cis  CH3ONO  -* 

CH30  +  NO 

41  .1 

53.1 

37.1 

33.3 

28.2 

Trans  CH^ONO 

CH30  +  NO 

41  .1 

38.1 

33.0 

27  .3 

Table  I.  Comparison  of  computed  activation  energies  and  experimental 
values.  Activation  energies  are  in  kcal/mole. 


At  (k  c  a  I ,  mole  J 


CH3N°2 


ch3  +  no2 


Experimental  Activation  Barrier 


B  WE  N 
B  WE  N 1 
B  WE  N2 
M  N  DO/CI 


C  N  BOND  L  E  NGTH ( A  ) 

Figure  3.  Computed  reaction  path  curves  for  the  dissociation 
of  n i tromethane .  The  arrow  indicates  the  vertical 
position  of  the  experimental  activation  energy. 


Trans  CHgONO 


7  0 


6  0 


5  0 


CH30  +  NO 
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01 

o  4  0 
E 


<0 

o 
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UJ 

<  30 


1  0 


1.0 


ON  BOND  LENGTH(A) 

Figure  2.  Computed  reaction  path  curves  for  the  dissociation  of 
trans  methyl  nitrite.  The  arrow  indicates  the  vertical 
position  of  the  experimental  activation  energy. 
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compared  these  with  Engelke's  MNDO/CI  results  and  with  experimental  vaiies. 
The  MNDO/CI  procedure  used  by  Engel ke  is  MNDO  plus  a  minimal  Cl,  where  tin’s 
minimal  Cl  includes  only  the  lowest  double  excitation.  The  calculated  BWEN, 
BWEN1 ,  and  BWEN2  reaction  path  curves  for  the  dissociation  of  CH^ONO  (cis  and 
trans)  and  CH^NO^  are  shown  in  Figs.  1-3,  respectively.  The  MNDO/CI  curves 
we  calculated  are  shown  for  purposes  of  comparison.  They  are  similar  to  the 
results  Engel ke  obtained.  Table  I  gives  our  computed  BWEN,  RWEN1  ,  and  BWEN2 

activation  energies  and  compares  them  with  Engel ke's  MNDO/CI  results  and  with 

20 

experiment  .  As  shown  by  the  table,  the  BWEN  activation  energy  for  CH^NC^ 
is  within  11  kcal/mole  of  the  experimental  value  and  is  only  a  slight 
improvement  over  the  MNDO/CI  result.  For  CH^ONO,  however,  the  BWEN  calculation 
is  a  substantial  improvement  over  the  MNDO/CI  result,  giving  activation 
energies  that  are  within  4  kcal/mole  and  3  kcal/mole  of  experiment  for  the 
cis  and  trans  stereoisomers,  respectively.  Engel  ke  calculated  the  MNDO/CI 
energy  only  for  cis  CH^ONO.  For  both  CII^ONO  and  CH^NO^,,  the  BWEN1  activation 
energies  are  not  as  close  to  experimental  values  as  the  BWEN  results  are, 
while  the  BWEN2  activation  energies  are  even  farther  from  experiment  than  the 
BWLN1  results.  Thus  for  the  activation  energies  of  methyl  nitrite  and 
nitromethine,  BWEN  gives  the  best  results.  In  both  cases,  these  are  improve¬ 
ments  over  the  MNDO/CI  results. 

Since  BWEN  gives  better  results  than  BWEN1  and  BWEN2  for  the  activation 
energies  of  CH^ONO  and  CH^NO^,  it  appears  that  the  procedure  that  works  best 
f r  larger  molecules  is  standard  second-order  perturbation  theory  in  which 
trie  zero-order  reference  state  is  the  HF  ground  state.  When  a  different  zero- 
order  reference  state  (the  lowest  root  of  a  two-configuration  Cl)  is  used  but 
the  excitations  are  still  with  respect  to  the  HF  ground  state,  as  in  BWENl 
and  BWFN2,  the  results  do  not  appear  to  be  as  good.  This  may  be  due  to  the 
inconsistency  of  using  excitations  from  the  HF  ground  state  while  using  a 
zero-order  reference  state  that  is  not  the  HF  ground  state.  If  this 
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over  o.l  states  that  are  doubly  excited  with  respect  to  the  HF  ground  state. 
The  MNDOC  procedure  also  contains  two  variants  of  BWEN,  called  BWENl  and 
BWEN2.  In  the  BWENl  method,  the  zero-order  reference  state  for  the  second- 
order  perturbation  treatment  is  not  the  ground  state  but  instead  is  the 
lowest  root  of  a  two-configuration  Cl  matrix,  where  the  two  configurations 
are  the  (IF  ground  state  and  the  lowest  doubly-excited  state.  The  perturbation 
expansion  then  extends  over  all  configurations  that  are  doubly-excited  with 
respect  to  the  HF  ground  state.  In  the  BWEN2  method,  the  zero-order  reference 
state  is  again  the  lowest  root  of  the  two-configuration  Cl  matrix  mentioned 
above,  but  the  perturbation  expansion  extends  over  single,  double,  triple  and 
quadruple  excitations  with  respect  to  the  HF  ground  state.  Since  the  BWEN2 
perturbation  expansion  extends  over  a  larger  configuration  space  than  BWENl 
does,  BWEN2  is  more  time-consuming. 

We  considered  the  following  reactions  for  the  dissociation  of  CH3N02 
and  CH3ONO: 

C!l3N02  -  CH3  +  N02 

CH3ONO  -  CH30  +  NO 

In  these  reactions,  the  dissociation  releases  highly  reactive  radicals  that 
then  attack  whatever  is  available,  releasing  more  energy  than  the  bond 
breaking  (or  dissociation)  required  and  thereby  producing  an  explosion.  In 
the  case  of  methyl  nitrite,  we  calculated  the  activation  energy  for 
dissociations  proceeding  from  both  the  cis  and  trans  configurations.  As 
methyl  nitrite  probably  contains  both  configurations,  the  dissociations 
could  proceed  either  from  both  of  them  or  from  just  one,  with  the  other 
f irst  converting  to  it. 

For  the  three  reactions  shown  above,  we  calculated  the  activation 
energies  of  the  molecules  using  the  BWEN,  BWENl  and  BWEN2  procedures  and 
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III.  MflDOC  Procedure 

The  MNDOc"^  method,  which  is  a  new  correlated  version  of  the  MNDO^' 

method,  is  a  semi-empirical  procedure  for  including  correlation  in  calculations 

of  molecular  energies.  One  of  the  advantages  of  semi-empirical  methods  is 

that  they  are  very  fast.  Due  to  its  speed,  the  MNDOC  method  can  be  used  to 

calculate  potential  surfaces,  transition  states,  and  activation  energies  of 

large  explosive  molecules.  We  have  recently  used  the  MNDOC  method  to  calculate 
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the  activation  energy  of  methyl  nitrite  (CH^ONO)  and  nitromethane  (CH^NC^). 

In  the  case  of  methyl  nitrite,  the  computed  results  were  in  good  agreement 
with  experiment,  indicating  that  this  fast  method  has  the  potential  to  yield 
accurate  molecular  activation  energies.  The  nitromethane  studies  have 
suggested  a  modification  of  the  MNDOC  procedure  which  is  expected  to 
significantly  increase  its  accuracy  and  thereby  make  it  applicable  to  a  wide 
range  of  explosive  molecules.  We  have  already  begun  work  on  this  new 
modification,  which  is  the  incorporation  of  a  two-configuration  SCF  into  the 
MNDOC  procedure.  The  two-configuration  SCF  is  expected  to  work  well  for 
molecules  like  nitromethane  in  which  there  is  a  small  energy  separation 
between  the  ground  state  and  first  excited  state.  The  two-configuration  SCF 
is  also  expected  to  work  well  for  larger  explosive  molecules  that  contain 
the  nitro  group,  such  as  nitrobenzene  and  trinitrotoluene  (TNT).  The  MNDOC 
method  and  results  are  described  in  more  detail  below. 

Second-order  perturbation  theory,  which  is  very  fast,  forms  the  basis  of  the 
MNNOC  method.  The  MNDOC  procedure  uses  Br i 1 1 ou in-Uigner  perturbation  theory  with 
f.ps tci n-Nesbet  energy  denominators  (BWEN).  In  the  BWFN  method,  the  second- 
order  perturbation  theory  is  based  on  a  single  reference  determinant:  the 
Hurtree-Fock  (IIF)  ground  state,  and  the  perturbation  expansion  extends  over 
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WRpA(!2)  =  W0(12)  .  (25) 

As  shown  by  Eqs.  (17),  (18),  (20),  (21)  and  (24),  the  wave  functions  and 
energies  that  appear  in  the  RPA  formula  for  the  dielectric  function  are 
Hartree  wave  functions  and  energies.  Therefore  one  should  use  Hartree  wave 
functions  and  energies  to  compute  RPA  dielectric  functions. 

Continuing  the  iterations  with  step  5),  we  put  I  (J  arid  inf  .,  iq.  (1) 
to  obtain  the  first  order  self-energy  and  first  or  for  Green's  function  G,  : 

(42)  =  ijW0(l+3)G1(14)rC|(42,';)d(.iT)  .  (26) 

’Jsing  [q.  (13),  Eq.  (22)  becomes 

2-j  ( 1  2 )  =  iG1  (12)WQ(1  +2)  .  (27) 

This  expression  for  is  the  first  term  in  the  expansion  of  in  powers  of 
the  screened  interaction  W.  The  main  point  here  is  that  we  are  expanding  in 
a  screened  interaction  rather  than  a  bare  Coulomb  interaction.  The  screened 
interaction  makes  it  possible  to  include  more  correlation  while  using  fewer 
terms  in  the  expansion  (compared  to  the  usual  type  of  bare  Coulomb  interaction 
expans  ion) . 

If  we  want  to  work  in  the  one-electron  picture,  we  must  obtain  an  expression 
fur  :: |  that  is  independent  of  frequency.  We  will  then  put  this  frequency- 
independent  :,i  into  Eq.  (13)  to  obtain  the  one-electron  wave  functions  u^  •  and 
energies  !\  .  that  are  used  to  define  (according  to  Eq.  (12)).  In  order 
to  obtain  a  self-energy  5.^  that  is  independent  of  frequency,  we  write  the 

p 

screened  interaction  Wq  as  a  bare  Coulomb  part  plus  a  part  Wq  due  to 
polarization : 


W0(l+2)  =  v(  1  +2)  +  W  q  ( 1  +2) 


(28) 


Thereto  re 


It  can  be  shown  that  1^(12)  is  a  sharply  peaked  function  of  t-j^. 

we  will  approximate  in  Eq.  (28)  by  its  integrated  value  times  a  I  function; 

however,  instead  of  6( ti -t^  +  n)  we  will  use  dtt-j^)  so  as  to  pick  up 

21  -23 

contributions  from  both  parts  of  the  Green's  fiction  ; 

w[j(l+2)  =  wj(x1,x2,u=0)6(trt2)  .  (29) 

Z-j  then  becomes  the  screened  exchange  plus  Coulomb  hole  (SECH)  expression  for 
21  -23 

the  self-energy  ,  which  contains  both  exchange  and  correlation: 


ESECH^X,X')  =  -P1  (x,x  ’  )WQ(x ,x  1  ,;u=°) 


6(x,x 1 ) [WQ(x,x 


--0) - . ; r, r  '  )  ] 


(30) 


where 


W0(x,x',u=0)  =  v(r.r') 


p , 

+  «  q  v  x ,  x  ,  .  - '  ! 


(31  ) 


and  p-j  is  the  density  matrix  computed  with  the  wave  functions  u^  .  : 

P1  (x,x  ' )  =  j:  u1  .  (x  1  ) u -j  •  (x)  .  (32) 

i  occ. 

According  to  Eq.  (22),  we  can  write  Wj  in  terms  of 


^0^xi  /  J  ( r  i  ,  r  ^ q  (x^iX^i-j'OjdXj 

wnere  cQ(x^  ,X2,a,=0)  is  the  static  RPA  dielectric  function.  As  noted  ear"! 

the  RPA  dielectric  function  corresponds  to  the  Hartree  equation  and  is 

calculated  with  fiartree  wave  functions  and  energies. 

Putting  which  is  independent  of  frequency,  into  Eq.  (13),  we 

o  L  L  •  1 

obtain  the  SECH  one-electron  equation: 


33' 


ler , 


23 


1 


[h(x)+V(x)]u1  .(x)  +  /xs[Cl|(x,x' )un(x' )dx'  =  F^j^lx) 


which  is  solved  sel f-cons i stentl y  to  obtain  correlated  wave  functions  and 
energies.  During  this  self-consistent  process,  the  screening  function  cq 
remains  constant  (it  is  not  recalculated  each  time  with  SECH  wave  functions  and 
energies,  but  rather  it  is  calculated  once  and  for  all  with  Hartree  wave  functions 
and  energies).  The  SECH  wave  functions  u^  •  and  SECH  energies  E-j  .  are  used  to 
define  the  first  order  Green's  function  G-j  ,  according  to  Eq.  (12): 


G,  ( x , x  1 ,  w)  =  l 


^x'  )u1  j(x) 


« 


r  ■  *  ’V" 

'.i 

C.  SECH2  fc|U3  t  i u n 

i 

j 

« 

f 

Putting  and  G,  into  Eq.  (2),  we  obtain  an  equation  for  the  first 

PP- 

order  vertex  function  I’-, : 

9 

6E,(12) 

I’,  (1  2,3)  =  6(12)6(13)  +  Gi(46)i'1(67,3)G1(75)d(4567  ). 

(36)  i 

l 

Then  using  Eq.  (27)  for  I,,  Eq.  (36)  becomes 

* 

r1  (12,3)  =  6(12)6(13)  +  iWQ(  1  +  2) J G-j  ( 1  6)  i'i  ( 67 ,3)G-j  ( 72)d( 67 ) . 

(37) 

Solving  the  above  equation  for  ,  we  obtain  an  infinite  power  series  in 

W0: 

► 

« 

1 1  (12,3)  =  6(12)6(13)  +  iW0(l+2)G1(l3)G1(32) 

-  W0(l+2)/G1(63)G1(37)WQ(6+7)G1(l6)G1(72)d(67) 

c 

-  iW0(l+2)/G1(83)G1(3y)W0(8  +  9)G1(tt)G1(97)WQ(6+7)G1(l6)G1(72)d(  6789)  j 

+  .... 

(3P) 

Tr, ;  s  and  explicitly  satisfy  the  Ward  identity.  Putting  Eq.  (38) 

;  into  Eq.  (3),  we  obtain  the  first  order  polarization  P, : 

for  j 

> 

m 

8,(12)  -  -10,(21)0,(12) 

i 

f/G, (bl )G1(16)W0(6f6)G,(2S)G,(62)d(56) 

- 

4 

+ i / G , (71 )G,(18)W0(7+8)G,(57)G1(86)W0(5+6)G,(26)G1(62)d(5678) 

I 

► 

-  JG, ( 91  )  G,  ( 1  , 1 0 ) WQ ( 9+l  0 )  G,  ( /  9 ) G ,  ( 1  0 , 8 ) W 0 (  7  * 8 ) G,  ( 6 7  ) G,  ( 86 ) WQ ( 5+ 6 )  - 

% 

x  G,  ( 25 )G,  ( 62  )d ( 567891  0) 

1  1 

\ 

f  .... 

(39)  : 

« 

i — *•-  - 

> 

.  ..... .  .  ... 

[1] 


The  above  expression  for  Pp  like  Eq.  (38)  for  ip  is  an  infinite  power 
series  in  Wg.  Putting  P-j  into  Eq.  (4),  we  obtain  the  corresponding  first 
order  dielectric  function  Ep 


£■,(12)  =  6(12)  -  /P1(13)v(32)d(3)  . 

(40) 

Then  putting  into  Eq.  (5),  we  obtain  the  first  order 

Up 

W1 (1 2)  =  /v( 1 3) e”1 ( 32)d( 3)  . 

screened  interaction 

(41) 

Since  the  above  expression  for  P-|  goes  beyond  the  RPA  to  include  higher  order 

correlation  terms,  we  call  it  the  generalized  RPA  (GRPA)  polarization,  and  we 

call  c i  and  W,  the  GRPA  dielectric  function  and  GRPA  screened  interaction 

respectively: 

W2)  '  pd'2> 

(42) 

"GP.PA^12^  =  C1  (12) 

(43) 

WGRPA^  ^  2^ 

(44) 

According  to  Eqs.  (39)  and  (42),  the  GRPA  formula  for  the  polarization  contains 
both  Wg  and  Gp  As  mentioned  previously,  Wg  is  the  RPA  screened  interaction 
and  is  computed  with  Hartree  wave  functions  and  energies,  while  G^  is  computed 
with  SFCH  wave  functions  and  energies,  according  to  Eq.  (35).  Thus  both 
Hartree  wave  functions  and  energies,  and  S E CM  wave  functions  and  energies  are 


Eq.  (39)  can  be  written  as  an  integral  equation  if  one  defines  the 


two -body  polarization  tt  by 


*0,2, 3, 4)  =  ttq(  1  ,2,3,4) 

./n(5,2,3,6)W0(5+6),0(l,5,6,4)d(S6) 

where 

nQ(l ,2,3,4)  =  -i61(12)G1(34)  . 

Solving  Eq.  (45)  for  v  yields 


t:C1  .2,3,4)  =  tt(j(  1  ,2,3,4) 

-/tt0(5,2,3,6)W0(5+6)t.0(1  ,5,6,4  )d(  56) 

+  /H0(7,2,3,8)Wc(7+8)’i0(5>7,d,6)W0(5  +  C)1>0O  ,5, 
-f-0(9,2,3JO)W0(9+10).0(719,10<P,)W0(7+8)TT0(£ 

l 

X  Wq(  5  6)  TTj, 


; f  o no  de  t  i nes 


P 1 ( 1 2 )  ;  *  (2,1  ,1  ,2) 


.mi  uses  •  j.  (47),  one  obtains  the  following  equation  lor 


(45) 

(46) 

6,4 )d( 5678 ) 

,7,8,6) 

,(1 , 5 ,6,4  )d(  567891  0) 

(47) 


"0uj  j  :n 


P,  (12)  = 

- J"q(5,1  ,1  ,6)W0(rjfG)  ^(^/j.G.L^dCjG) 

+  1^(7  ,1  ,1  ,8)Wq(  7+8 )  ttq(  5, 7,8,6)Wq(  5  +  6)  7rQ(2,b,6 ,2)d(  5678) 

- /ttq(9,1  ,1  ,1  0)W0(9+10)ttq(7,9,10,8)W0(7+8)tt0(5,7,8,6) 

x  W0(5+6)TTQ(2,5,6,2)d(5678  91  0) 

+  ....  (49) 

Then  putting  Eq .  (46)  for  Uq  into  Eq.  (49),  one  obtains  Eq.  (39).  Thus  (1 2 ) 
is  given  by  the  integral  equation  (45)  if  one  first  solves  Eq.  (45)  for  v 
and  then  uses  Eq.  (48). 

One  can  solve  Eq.  (45)  for  tt  in  terms  of  inverses  if  one  defines  the 
inverses  -  \  ^  and  by 


/:(1  ,2,3,4)r1(4,6,5,l)d(14)  = 

■5(25)6(36) 

(50) 

/  rr  ( 1  ,2,3,4b'1  ( 6,3 ,2,5)d(23)  = 

6(15)6(46) 

(51) 

/■i(l  ,2,3,4)  ~  ~ 1  (4  ,6,2 , 5 )  d  ( 24 )  = 

6(15)6(36) 

(52) 

/t:(1  ,2,3.4)tt-1  (6,3,5.1  )d(13)  = 

6(25)6(46) 

(53) 

--1 (1  ,2,3,4)  =  iGj1 (12)G^ (34) 

(54) 

/G1(12)G[1(23)d(2)  -  £(13) 

(55) 

fG,  (12JG71  (31  )d(l  )  =  5(23) 

(56) 

where,  according  to  Eqs.  (46)  and  (54)-(5G), 

/~0(1  ,3.3, 4),’1(4, 6,5,1  )d(14)  =  6(25)6(36) 


(57) 


/iiQ  ( 1  ,2,3,4)irj^(6,3,2,5)d(23)  =6(15)6(46)  ((j8) 

Jtt^I  ,2,3, 4)*"1  (4,6,2, 5)d(24)  =  6(15)6(36)  '  (09) 

/1r0(l,2,3,4)1i',(6,3,5,l)d(  1  3)  =  6(25)6(46)  .  (  60) 

Multiplying  both  sides  of  Eq.  (45)  by 
,t'1(10,3,2,7)tt'1(4,9,8,1)  , 

integrating  over  coordinates  1,2,3  and  4,  and  then  replacing  coordinates 
10,9,8,7  by  1,2, 3, 4  respectively,  one  obtains 

-t'1  ( 1  ,2,3,4)  T,^7  (1  ,2,3,4)  +  W0(3  +  l  )o(  12)5(34  )  .  (61) 

Using  Eq.  (54)  for  it"1  ,  Eq.  (61)  becomes 

-r^O  ,2,3,4)  =  6^(12)G^(34)  +  WQ( 3  +  l  ) 6(  1  2) 3(  ,:4 )  .  (62) 

is  then  given  by 

■' ( ’  ,3,3,4)  =  (1  ,2,3,4)3“1 

=  [G^(12)G'1(34)  +  VJ0(3  +  1)6(12)6(34)]"1  (63) 

■vhere  the  inverse  is  defined  by  Eqs.  (50)  -  (53).  After  t ( 1  ,2,3,4)  has  been 
calculated  from  Eq.  (63),  P^(12)  is  then  given  by  Eq.  (48). 

Continuing  the  iterations  with  step  9),  we  put  !]  and  W]  into  Eq.  (1)  to 
determine  ;  and  G2 : 


34 


4 

f 

51 

I 

i 

€ 

i 


4 


E2(12)  =  i/W1(1+3)G2(14)r1(42,3)d(  34  ).  (  64) 

Putting  Eq.  (38)  for  r -j  into  Eq.  (64),  we  obtain  an  infinite  series  for 
::2(1  2)  =  iG2(12)w1(l+2) 

-/G204)W1(l+3)G1(43)G1(32)WQ(4+2)d(34) 

-i /G2(14)W](1+3)G1(63)G1(37)W0(6+7)G1 (46)G1(72)W0(4+2)d(3467) 

+  /G2(14)W1(1+3)G1(83)G1(39)W0(8+9)G1  (b8)G1(97)W0(b+7)G1(4  6) 

x  G] ( 72)W0(4+2)d (346789) 

+  .  ..  (65) 

,,'e  now  define  to  be  the  first  two  terms  in  the  above  series: 

::^(12)  =  iG^(l  2)W1  (1  +2) 

-/G2(14)W1(l+3)G1(43)G1(32)W0(4+2)d(34)  .  (66) 

The  above  expression  for  l 2  is  the  first  two  terms  in  the  expansion  of  l 

in  powers  of  the  screened  interaction  Vi,  and  62  is  the  corresponding  Green's 

function.  We  will  call  b2  the  second  order  self-energy  since  it  contains  the 

2 

first  and  second  terms  of  the  power  series  (it  contains  the  W  term  and  the  W 
term),  and  we  will  call  G2  the  second  order  Green's  function. 

In  order  to  continue  working  in  the  one-electron  picture,  we  must  obtain 
an  expression  for  Z2  that  is  independent  of  frequency.  We  will  then  put  this 
into  Eq.  (13)  to  obtain  a  new  one-electron  equation  that  contains  higher  order 
correlation  terms.  This  new  one-el ec Iron  equation  will  yield  correlated  wave 
functions  u2l-  and  correlated  energies  E2-  which,  through  Eq.  (12),  will  yield 


« 


the  second  order  Green's  function  G^.  In  order  to  obtain  a  frequency- 
independent  expression  for  we  treat  the  screened  interactions  that  appear 
in  Wg  anc*  wl  )  the  same  manner  as  Wq  was  treated  in  the  derivation  of 

the  SECii  equation.  In  other  words,  Wq  and  W^  are  written,  according  to 
Eq.  (28),  as  a  bare  Coulomb  part  plus  a  part  due  to  polarization: 


Wq(4+2)  =  v(4+2)  +  wJ(4+2)  (67) 

W-,  (1 +3)  =  v(  1  +3)  +  wfO+3)  .  (68) 

P  P 

Then  we  approximate  Wq  and  W,  by  their  integrated  value  times  a  6  function: 

Wq(  4  +  2 )  =  WQ(x4,x2>w=0)i(t4-t2  hi)  (69) 

wf(l+3)  =  w5>(x1,x3,w=0)6(t1-t3  +  rl)  .  (70) 


In  some  cases,  we  remove  the  n  from  the  6  functions  in  Eqs.  (69)  and  (70) 
so  as  to  pick  up  contributions  from  both  parts  of  the  Green's  function, 
according  to  Eq.  (29).  In  other  cases,  we  retain  the  n  in  the  i  functions  in 
Eqs.  (63)  and  (70)  so  as  to  pick  up  contributions  from  only  one  of  the  parts 
of  the  Green's  function.  Since  each  of  the  Green's  functions  in  Eq.  (66) 
has  two  parts  (c.f.  Eq.  (8)),  then  the  second  term  of  Eq.  (66)  yields  eight 
terms  when  Eq.  (8)  is  substituted  for  G.  By  carefully  considering  the  time 
dependence  of  the  Green's  functions  in  each  of  these  terms  and  by  using  the 
appropriate  functions  in  each  case,  we  obtain  the  following  frequency- 
independent  expression  for  which  we  call 


i'SECH2^xl  ’ x2^  =  -p2^xl  ,x2^1  ^X1  >x2,w-^ 

+  ~2  o( x-j  , x 2 )  [W ^  (Xi  .x^.w-O)  -  v( >  ^ ,  r2 )  ] 

- 2 ^ x i  >x4 )^i  ( X-j  ,  x^ , ) P]  (x4«x3)p-)  (X3»X2^0^X4’X2’ ) dx^dx4 
+  2  /p2^xl  ,x4 ^xl  ,x3* uj-^ ^ pi  (x4 >x3 )  x3 > x2 )  ^q(  x4  ,  x2 , w-  0)-v(r'(^,f'2)^x3^x4 
+  2'  / 5 ( x i  ,x^  )  [W-|  (xi  ,x3, w=0 ) -  V ( r-j , r^ ) 3 P]  ( x4  •  x3 ) P]  ^ x3 * x2 ^0 ^ x4  ’ x2 ’ ) d x 3d X4 
-  4  /fi(xpx4  )  [Wj  (x-j  ,X3 ,-x'  r  0) - v ( r]  »r3 ) ] P 1  (X4  ^3)  <$( .x^ )  [Wq(x4  ,x2,u>=0) 

-v(rr4,r2)]dx3dx4 

where  p2  is  the  density  matrix  computed  with  the  wave  functions  u^.  : 

p 2 ( x , x 1 )  -  z  u2i (x ' )u2i (x)  •  (72) 

i  occ . 

As  shown  by  Eq.  (71),  contains  four  higher  order  correlation  terms  (in 

addition  to  the  screened  exchange  and  Coulomb  hole  terms).  The  first  of  these 

higher  order  terms  is  the  product  of  two  screened  exchanges,  the  second  is  the 
product  of  a  screened  exchange  and  a  Coulomb  hole,  the  third  is  the  product  of  a 
Coulomb  hole  and  a  screened  exchange,  and  the  fourth  is  the  product  of  two 
Coulomb  holes.  Thus  "-ojrQ^  is  a  natural  second  order  extension  of  the  screened 
exchange  plus  Coulomb  hole  self-energy  v 

Putting  -5 2 CH2  into  (^),  we  obtain  a  new  one-electron  equation  that 
goes  beyond  the  SCCH  equation  to  include  higher  order  correlation  terms: 

[  h  ( X  )  f  V  ( x )  ]  a  i  ( x  )  4  /  f.s  E  CH  2  ( x , x  •  )  U  2  i  ( x '  ) d  x  1  =  E2iu2i^x^  •  (73) 

We  call  this  new  equation  the  SECH2  equation.  It  is  solved  sel f-consistently 
to  obtain  wave  functions  u2-  and  energies  E2^  that  have  a  high  degree  of 
correlation.  During  this  self-consistent  process,  the  screening  functions 


37 


lq  and  c -j  remain  constant  (they  are  not  recalculated  cacti  time  with  SECH2 
wave  functions  and  energies,  but  rattier  is  calculated  once  and  for  all 
with  Marti  ee  wave  functions  and  energies,  and  ^  is  calculated  once  and  for 
all  with  SF.CH  wave  functions  and  energies).  Also  p-j  remains  constant  during 
the  self-consistent  process  (p-j  is  calculated  once  and  for  all  with  SL'CH  wave 
functions).  The  only  quantities  that  change  during  the  self-consistent  process 
are  the  SECH2  wave  functions  u^ •  and  the  SECH2  energies  E^  • . 

Ttie  SECH2  equation  is  a  new  one-electron  equation  that  has  never  been 
used  before  in  either  molecular  or  solid  state  calculations.  It  is  expected 
to  lead  to  substantial  improvements  in  the  calculation  of  correlated 
molecular  energies. 
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